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We introduce the two-dimensional Gross-Pitaevskii/nonlinear-Schrodinger (GP/NLS) equation 
with the self-focusing nonlinearity confined to two identical circles, separated or overlapped. The 
model can be realized in terms of Bose-Einstein condensates (BECs) and photonic-crystal fibers. Fol- 
lowing the recent analysis of the spontaneous symmetry breaking (SSB) of localized modes trapped 
in ID and 2D double- well nonlinear potentials (also known as pseudopotentials), we aim to find 
2D solitons in the two-circle setting, using numerical methods and the variational approximation 
| (VA). Well-separated circles support stable symmetric and antisymmetric solitons. The decrease 

, of separation L between the circles leads to destabilization of the solitons. The symmetric modes 

undergo two SSB transitions. First, they are transformed into weakly asymmetric breathers, which 
is followed by a transition to single-peak modes trapped in one circle. The antisymmetric solitons 
perform a direct transition to the single-peak mode. The symmetric solitons are described reason- 
ably well by the VA. For touching (L = 0) and overlapping (L < 0) circles, single-peak solitons 
are found — asymmetric ones, trapped in either circle, and symmetric solitons centered at the mid- 
point of the bi-circle configuration. If the overlap is weak, the symmetric soliton is unstable. It 
may spontaneously leap into either circle and perform shuttle motion in it. A region of stability of 
\yy ' ' the symmetric solitons appears with the increase of the overlap degree. In the case of a moderately 
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strong overlap, another SSB effect is found, in the form of a pair of symmetry-breaking and restoring 
bifurcations which link families of the symmetric and asymmetric solitons. 

Keywords: soliton; breather; Bose-Einstein condensate; photonic-crystal fiber; Gross-Pitaevskii 
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I. INTRODUCTION 

> \ 

A large number of systems in optics and atomic physics feature a dual-core (double-well) structure, based on an 
effective potential in the form of a symmetric set of two minima (wells). The wells are linked by the tunneling of 
optical waves or matter waves, trapped in them, across the barrier separating the wells. In cases when the trapped 
. ■ waves are subject to nonlinear self-interaction, a fundamental property of the symmetric double-well systems, which 
originates from the competition between the linear tunnel coupling between the wells and the nonlinearity acting 
inside of them, is the spontaneous symmetry breaking (SSB). In particular, the SSB of solitons (solitary waves trapped 
in the double- well potential) was first studied in the model of dual-core optical fibers ill - [6l| . A similar analysis was 
carried out for Bose-Einstein condensates (BECs) loaded into double- well potentials [7|-[l0j. In the experiment, the 
SSB has been observed in the BEC with repulsive interactions between atoms [Tlj |. an d in an optical setting based 
on photorefractive crystals [12J. It has been demonstrated that, in the systems with attractive or repulsive intrinsic 
5_j ■ nonlinearities, stationary asymmetric modes are generated by symmetry-breaking bifurcations from symmetric or 
antisymmetric trapped solitons, respectively. SSB effects for solitons were also studied in the model with competing 
nonlinearities of both these typ es, viz., self- focusing cubic and self-defocusing quintic nonlinear terms, which gives rise 
to closed bifurcation loops [lj, [lj| . 

The analysis of the SSB in matter-wave (BEC) settings and their photonic counterparts was extended for two- 



dimensional (2D) solitons, supported by the interplay of the self-attractive [lE|-[T7l or self-repulsive [16|-[18| nonlin- 
earity and spatially periodic potentials acting in each core. The periodic potential is necessary to stabilize 2D solitons 
(in each core separately) against the collapse in the case of the self-attraction [l9[ , or support solitons of the gap type 
in the case of the self- repulsion [2(|-[22j]. The analysis of the SSB in a 2D localized configuration of a different type, 
based on the set of four potential wells forming a square pattern (in a single core) was reported in Ref. (23[. 

Photonic and matter waves can be trapped not only in the usual linear potentials, but also in nonlinear pseudopo- 
tentials (this term is often used in solid-state physics [24jh which may be induced by a spatial modulation of the local 
nonlinearity coefficient. In BEC, the local nonlinearity may be controlled, through the Feshbach-resonance effect, by 
an external magnetic field (25j . Accordingly, the spatial nonlinearity modulation may be induced by an inhomogeneous 
magnetic field [26]. In optics, similar settings can be designed as all-solid or liquid- filled microstructured fibers, using 
combinations of materials with matched refractive indices but different Kerr coefficients, as proposed in Ref. [271 ] . 

The studies of soliton dynamics in purely nonlinear (pseudo)potentials, as well as in their combinations with usual 
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linear potentials, has recently grown into a vast research area, see the review arti cle 12811 . Models with various patterns 
of the spatial nonlinearity modulation were studied in detail in ID geometries |29l|-[35j. [28j . In 2D, the analysis is 
essentially more challenging, as, in this case, it is difficult to stabilize localized modes against the spatiotemporal 
collapse (catastrophic self- focusing induced by the cubic nonlinearity) , using only the local modulation of the nonlin- 
earity, without the support provided by a linear potential [36j|-[38|], [27|], (28|. Nevertheless, 2D nonlinear structures 
capable to support stable solitons have been found. Basic types of such structures are represented by a circle or 
annulus, filled with the self-focusing material, which is embedded into a linear or defocusing medium [36j j (as well 
as a lattice of such circles [27[), and a single or double stripe made of the same material [38). In both cases, it was 
concluded that the nonlinearity-modulation profiles with sharp edges provide for much more efficient stabilization of 
2D solitons than similar smooth profiles. On the other hand, no pattern of 2D modulations of the self-focusing cubic 
nonlinearity was found that would support stable vortex solitons. 

Coming back to ID, it was found that the simplest pseudopotential capable to sustain stable localized modes is 
built as a symmetric set of two delta- functions multiplying the cubic self- focusing nonlinearity [33[ . Explicit analytical 
solutions are available for all stationary modes trapped in this setting - symmetric, antisymmetric, and asymmetric 
ones. The analytical solution allows one to study the corresponding SSB bifurcation, which generates asymmetric 
modes from the symmetric ones. On the other hand, the exact solution of the ID model with two delta-functions is 
partly degenerate, which, in particular, makes the asymmetric modes completely unstable. The replacement of the 
ideal delta-functions by their regularized counterparts lifts the degeneracy, making the asymmetric modes partly or 
fully stable H. 

The model with a pair of ideal or regularized delta-functions may be used as a paradigm for the study of the SSB 
of localized modes supported by symmetric double- well pseudopotential structures. The objective of the present work 
is to introduce a natural generalization of this setting in the 2D geometry, with the nonlinearity concentrated in two 
separated or overlapping identical circles with sharp edges, see Fig. [T]and insets to Figs. [Hfa) and[T27a) below (ideal 
(5-functions are irrelevant in the 2D case, see Eqs. <jTUJ) , (|TT|) and the related text in the next section). This setting 
can be realized in BEC and nonlinear optics alike, by means of the same techniques which were outlined above, i.e., 
respectively, the use of the Feshbach resonance, controlled by a spatially nonuniform magnetic field, or a crystalline 
structure formed by a pair if intermingled materials with equal refractive indices but different Kerr coefficients. 

The paper is organized as follows. The model is introduced in detail in Section 2. Then, in Section 3 we develop 
the variational approximation (VA), with the objective to describe patterns with different symmetries supported by 
the symmetric pair of well-separated circles. Results of the numerical analysis are collected in Section 4. The results 
describe the transition from double-peak symmetric and antisymmetric modes to asymmetric ones, with the decrease of 
the separation between the circles, which implies the strengthening interaction between them. The symmetric solitons 
perform the symmetry-breaking transition via an intermediate mode in the form of weakly asymmetric breathers, while 
antisymmetric solitons directly jump into strongly asymmetric single-peak modes. For touching or partly overlapping 
circles, only single-peak solitons exist. These may be either symmetric ones, centered at the midpoint of the bi-circle 
configuration, or solitons spontaneously leaping into either circle and then featuring shuttle motion in that trap. 
Various manifestations of the SSB in the present system are summarized in Section 4, which concludes the paper. 



II. THE MODEL 



A. The equations 



The setting considered in this work is shown in Fig. [TJ where a is the radius of the two identical circles, and L is 
the separation between them, R = L + 2a being the distance between the centers of the circles, which are located at 
points {xo, Ho} = {± (L/2 + a) , 0}. Along with the configuration displayed in Fig. [TJ we will also consider the one 
with R < 2a (i.e., L < 0), which corresponds to partly overlapping circles [see top left insets to Figs. and IT2T s) 
below] . 

The scaled form of the underlying Gross-Pitaevskii/nonlinear-Schrodinger (GP/NLS) equation with the self- focusing 
(.9i > 0) nonlinearity confined to the two circles is 



(ipxx + Ipyy) ~ 9 (X, y) |V>| ^, 



(1) 



where the simplest form of the localization may be chosen as 



9(x,y) = 9i 



{x + R/2) 2 + i/ 



(2) 
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FIG. 1: (Color online) The geometry of the bi-circle system. 



In the case of the GP equation, t is normalized time, while in optics, t is actually the propagation distance in the 
bulk waveguide. In either case, x and y are the transverse coordinates. In fact, a different form of the 2D modulation 
function will be used for numerical calculations, see Eqs. (j!?Tl -(l2"3"l below. 

Stationary solutions to Eq. (JIJ, with chemical potential [i of the BEC (or propagation constant — [i in the optical 
waveguide) are sought for as ip (x, y) = e~ vt (a;, y), where real function ip (x, y) obeys equation 

Ll<j)+ {<j)xx + <Pyy) + gix^)^ = 0. (3) 

Equation ^ will be solved below by means of a numerical method, and the VA will be used too. 



B. The Lagrangian structure and variational ansatz 



Equation §3§ can be derived from the Lagrangian, L = J J Cdxdy, with density 

2C = + {^f + {<t>yf - (V% (x, y) 4 . 
Obviously, localized solutions have fi < 0. The energy corresponding to Eq. (flj is 



E = 



\g (x,y) |-0 (x, y) | 4 



dxdy. 



(4) 



(5) 



The application of the VA is possible in the case when a in Eq. @ is small, hence g (x, y) may be approximated 
(only for the purpose of the development of the VA) by a combination of delta-functions: 



g (x, y) = G [5 (x - R/2) 6 (y) + 6 (x + R/2) S(y)] , G = 7r gi a 2 



(G) 



where coefficient G is determined by the integral-balance condition, J J g (x, y) dxdy = 
G J f [S (x — R/2) S (y) + 5 (x + R/2) 5(y)]dxdy. The VA developed below is based on the use of the following 
ansatz, 



(j> (x, y) = exp 



2y 2 

R?l 2 



A exp 



2(x-R/2) 2 
R 2 P 



Bexp 



2(x + R/2Y 
R 2 l 2 



(T) 



where A, B and I are variational parameters, the latter one being the width of the ansatz measured in units of R/2. 
The ansatz is isotropic in the vicinity of each attractive center (the variational equations for an anisotropic ansatz 
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turn out to be very cumbersome). It is assumed that the asymmetry between the wave functions in vicinities of the 
two attractive centers may be accounted for by a difference in the amplitudes, A ^ B, while widths m and I are 
assumed to be the same near both centers. 

Comparing expression with Eqs. ([2]) demonstrates that approximation ([5]), based on the delta- functions, may 
be valid for sufficiently small a in Eq. ([2]) — namely, if the VA will eventually yield 



I > 2a/ R. 



(8) 



On the other hand, it is easy to check that Eq. (j3|) with the nonlinearity-modulation function represented by the 
sing le 2D delta-function, i.e., g (x, y) = G5(x)8(y) 1 does not give rise to any soliton solution (unlike its ID counterpart 
|33|). Indeed, taking the corresponding equation, 

p<j) + (<j) ) + G6 (x)S(y) 4> 3 = 0, (9) 

and setting, at first, p — 0, one can use the commonly known fact that the Poisson equation in the form of 

4>xx + 4> V y = C5(x)6{y), (10) 

with constant C, has the fundamental solution, 



(r) = (C/27r)ln(r/r ) ; 



(11) 



where r = y 1 x 1 + y 2 , and ro is an arbitrary positive constant. Temporarily introducing a small cutoff radius, p, 
formula (fTTj) suggests that a solution to Eq. ((9|) can be sought for as cf)(r) — cf>oln(r/ro) , which corresponds to 
C = —G(<p(r = p)) 3 = G<pQ (ln(r/p )) 3 = C p in Eq. (TTOj) . Then, the self-consistency condition ensuing from Eq. 
(fTTj) . 0o = C p /2ir, yields 0§ = (2-k/G) [In (ro/ p)]~ .This relation shows that the limit of p —> corresponds to (j>o — > 0, 
which means that solely the trivial solution, = 0, is possible if the cutoff is removed. 

A similar conclusion [the nonexistence of nontrivial solutions to Eq. (J9j)] can be obtained for p < 0, using an 
appropriate Hankel's function, instead of solution (|11[) . Thus, the use of approximation © in the framework of the 
VA is meaningful for the circles of a small but hnite radius a, which obeys condition (f8]). 



III. THE VARIATIONAL ANALYSIS 

The substitution of ansatz ([7]) into density (j4j and calculation of the integrals yields the following expression for 
the effective Lagrangian: 

L cS = - (tt/8) R 2 pl 2 [A 2 + B 2 + 2e- 1/l2 AB"j 
+ (it/2) \(A 2 + B 2 )+2(1- r 2 ) e-^AB 
-G(R/4f [(l + e - 8 / ;2 ) (A 4 + B 4 ) 
+Ae- 2 ' 12 (l + e- 4 /' 2 ) AB {A 2 + B 2 ) + ^e^ 1 A 2 B 2 \ , (12) 

the norm of the ansatz being 

A = y J (f) 2 (x,y)dxdy = tt{RI/2) 2 (A 2 + B 2 + 2e- 1/l * AB^j . (13) 
In particular, for the symmetric and antisymmetric solitons, with A = ±_B, Eqs. (|12[) and (|13[) give 



(R/2) 2 pi 2 [1 ± e 
-(G/8)R 2 A i \l + e-^ 2 ± 4e- 2 / /2 



i±(i-r 2 ) e- 1 ' 12 



1 + e 



-A/V 



6e 



-A/l 2 



(14) 



N = rrl 2 (R 2 /2) (l + e^ 1 /' 2 ) A 2 



(15) 
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Even after the simplifications, the Euler-Lagrange equations, dL^'/d (A 2 ) = dL^'/d (l 2 ) = 0, following from Eq. 
(1141) remain cumbersome. A tractable approximation can be developed if I 2 is small enough to neglect exponentially 
small factors e~ 4 ^ 1 and e~ 8 /' in comparison with 1 in expression (TT41 (in fact, we will then have I 2 < 1/2, see below). 
In this approximation, the variational equations for the symmetric solitons are 



i + (r 2 + i) 



-(2/i?rr 4 (r 2 -2) e 



(16) 



A 2 



1 - (I- 2 1) 



(R/2) V/ 2 ( 



l + e- 1 ' 1 '' 



G(R/2) 2 (l + 4e- 2 /' 2 ) 



(17) 



It follows from Eq. (fT6|) that the localization condition, /i < 0, entails I 2 < 1/2, as mentioned above. Then, both 
terms in the square brackets in Eq. (|17|) are positive, hence A 2 is positive too, as it should be. 

An explicit analysis is also possible for the limit case of very broad solitons, I 2 3> 1. In this case, the expansion of 
Lagrangian (|T2j) in powers of l~ 2 yields 



L e f[ 



- (vr/8) i?V 2 (A + B) 2 

+ (tt/2) (A + B) 2 — G (R 2 /8) (A + B) 4 



-2irAB + G (R 2 1 2) {A + B) 



(18) 



A straightforward consideration of the the Euler-Lagrange equations following from Eq. (fTS]) demonstrates that they 
give rise to no antisymmetric solutions, while the symmetric solution is found in the following form, valid at the lowest 
order in l~ 2 : 



A 2 = B 2 = tt/ (2GR 2 ) , N 
N = (n 2 /GR)(-»y 1/2 . 



[n/(2G)] 2 l 2 , 



(19) 
(20) 



In this approximation, Eq. (|19[) demonstrates that the amplitude of the symmetric soliton is constant, while its norm 
diverges ~ I 2 , and the entire soliton family is described by dependence N(n) given by Eq. (|2U)) . 



IV. NUMERICAL RESULTS 



A. The setting 



For the numerical analysis, we return from approximation ([6]) to the configuration displayed in Fig. [I] and define 
the nonlinearity as follows [cf. Eq. ||2J)]: 



9i V) = ex P - i r L/a) k + exp - (r R /a) k 



(21) 



92 (x,y) 



exp 



{tl/ a) 



exp 



(22) 



g{x, y) = (1/2) [gi (x, y) + g 2 (x, y)] , 



(23) 



where Tl = \j y 2 + (x + R/2) 2 and r R = \J y 2 -\- {x — R/2) 2 . As mentioned above, the stabilization of 2D solitons 

supported by the single nonlinear circle is facilitated by making edges of the circles sharp [27], [HI, [3(1 . To implement 
this feature in the present setting, k = 100 was used in Eqs. (l2~Tj) and (|2"2"|) [the full nonlinearity (|2"3"|) is then almost 
identical to g\ (x, y)]. 

Localized stationary modes were found as solutions to Eq. (j3]) by means of the so-called Newton-conjugate-gradient 
method [3{| (the ordinary version of the Newton's method did not lead to convergent solutions in the present setting) . 
After finding the stationary solutions, their stability was tested by means of direct simulations of Eq. (TT]), performed 
by means of the usual split-step method. 

In addition to presenting the numerical results, we also include some predictions produced by the VA [chiefly, in 
Fig. Ufa)]. To this end, norm (IT5|) was plotted vs. /i, taking A 2 as per Eq. ()17|) and eliminating I, in favor of /z, with 
the help of Eq. (E 
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B. Separated circles (L > 0) 

1. Symmetric solitons 

We start by presenting results for symmetric modes, with <j>(x,y) — <p(—x,y), in the case of L > 0, i.e., non- 
overlapping circles, see Fig. [TJ Figures (Ha) and (b) display the overall characteristics of families of the symmetric 
modes, in the form of dependences N (/i) and E(N), i.e., the norm versus the chemical potential and the energy versus 
the norm [see Eqs. (fT5|) and ©], for different values of radius a of the circles and a fixed separation between them, 
L = 3-7T (this value was selected to display generic results). 

The curve marked "a = 0" in Fig. [Ha) displays the prediction of the VA for the symmetric solitons, obtained from 
Eqs. (fT5|) . ([17]). and (|16p . as outlined above [recall that the VA was developed for the nonlinearity modulation taken 
in the form of Eq. ([6]), which formally corresponds to a = 0]. The shape of the VA-predicted curve is reasonably close 
to its numerical counterparts generated for finite a. The other tractable approximation, which amounts to Eq. (|20[) 
for [i — > —0 and N — > oo, correctly predicts the asymptotic form of the numerically found curves in this limit. 

The opposite limit of /i — » — oo corresponds to very narrow peaks pinned to the center of each circle. When the 
radius of these peaks is much smaller than the radius of the circles (a), they are nearly identical to the Townes solitons 
[4fj| . whose norm (i.e., the collapse threshold in the uniform medium with the self-focusing cubic nonlinearity) is 

N tht a 11.69. (24) 

This fact explains the limit value N(fi — > — oo) « 23.5 ~ 2A t hr observed in Fig. [UJa) for a = 1, 1.5, 2, and 3. The 
form of the VA adopted above is irrelevant in this case, as approximation ([6]) implies that the radius of the peak is 
much larger than a. However, the isolated peak in the medium with the uniform cubic nonlinearity (i.e., the unstable 
Townes' soliton) may be approximated by another version of the VA, which yields a reasonable approximation for the 

collapse threshold, N^ r A) = 4tt 

The plots displayed in Figs. [2fa) and[2jb) identify the stability of the families of symmetric solitons. A well-known 
necessary, but not sufficient, condition for the stability of solitons supported by the self- focusing nonlinearity amounts 
to the Vakhitov-Kolokolov (VK) criterion, dN/dfj, < (4(|. It is seen in Fig. HJa) that this criterion is indeed necessary 
but not sufficient, as the portion of the N(/j,) curve with the negative slope is stable only for a = 1. It is worthy to 
note that, as seen in Fig. [2] the stable portion of the latter curve corresponds to lower values of the energy than the 
unstable part of the same curve. 




H N 

(a) (b) 

FIG. 2: (Color online) Norm N versus fi (a) and energy E versus iV (b) for symmetric solitons at different values of radius a 
of the circles, with a fixed separation between them, L = 37r. The stability of the solution families is identified as follows: the 
magenta color pertains to the instability through decay; black - the instability via a jump into a strongly asymmetric mode; 
red - stability (it is a part of the curve for a = 1 with the negative slope); blue - spontaneous transformation of the stationary 
solitons into robust breathers. 

A typical example of a stable symmetric (double-peak) soliton is displayed in Fig. [3] (it pertains to the separation 
between the circles L = R — 2a = 6tt — 3 « 15.85, which is different from L = 3ir fixed in Fig. [5J but there is no 
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essential difference in the shape of the stable solitons in these cases). The evolution of weakly unstable symmetric 
solitons leads to the formation of breathers, as shown in Fig. |U In this case, the soliton is spontaneously transformed 
into a robust localized object featuring irregular small-amplitude intrinsic oscillations. It is worthy to note that the 
transition from the stationary solitons to breathers leads to a relatively weak but tangible SSB effect, as seen in Fig. 
Etb). Numerical data demonstrate that the breather keeps practically the entire initial norm, i.e., radiation losses are 
negligible in the course of the rearrangement of the symmetric soliton into the breather. 




FIG. 3: (Color online) An example of a stable symmetric soliton, with norm N = 22 (/i = —0.23045), R = 6ir, and a — 1.5. 
The evolution of the soliton is shown in the x-cross section drawn through y — 0. 

The instability mode designated as the decay into radiation in Fig. [2]means that the soliton completely disintegrates 
(not shown here in detail). On the other hand, the evolution outcome marked as the instability via a jump into a 
strongly asymmetric mode implies that the symmetric double-peak soliton suffers strong SSB, which leads to a sudden 
transition into a robust single-peak mode nested in either circle, as shown in Fig. [5] Note that, in the case displayed 
in Fig. El the initial norm (N — 22) is much higher than critical value ([24]) . The numerical data demonstrate that 
the eventual value of the norm of the emerging single-peak mode is slightly smaller than the threshold value (24), the 
surplus norm, AiV ~ 10, being shed off in the form of radiation. 

The results of the numerical analysis of the stability of the symmetric solitons is summarized in the diagram in the 
plane of (a,R/ir), which is displayed in Fig. [6] for a fixed total norm, N — 22. An obvious feature is that radius 
a of the circles supporting stable symmetric solitons must not be too small (roughly speaking, all the solitons are 
unstable at a < 1). In fact, stable symmetric solitons (as well as antisymmetric ones, dealt with in the next subsection, 
see Fig. [5] below) may be considered as pairs of virtually non-interacting stable fundamental solitons supported by 
isolated nonlinear circles - essentially the same fundamental modes which were reported in Refs. [36j and [27]. On 
the other hand, the further increase of a and/or decrease of R lead to the destabilization of the symmetric soliton 
- first, through its conversion into the breather (see Fig. @|, which is followed by the sudden transformation of the 
breather into a strongly asymmetric single-peak soliton (see Fig. [S]). These transitions can be readily understood: the 
increase of a and/or decrease of R lead to the enhancement of the interaction between the two circles, which causes 
the destabilization of the double-peak soliton. As concerns the stability border of the stable symmetric solitons, it 
has been checked that it is precisely determined by the VK criterion, i.e., it coincides with the locus of points where 
dN/dn — [recall, however, that these points do not always determine a stability border, as soliton families which 
include portions with dN/dfj, < may be unstable, see Fig. Eta)]. 

2. Antisymmetric solitons 

Figures [7] and |S] are counterparts of the above figures E] and [fj] for antisymmetric solitons, with <p (x, y) = —<f> (—x, y). 
Predictions of the VA for antisymmetric solitons are not included into Fig. Eta), as, on the contrary to the symmetric 
solitons, the agreement between the VA and the numerical findings is poor in this case. The diagram in Fig. 0includes 
a region of stable antisymmetric solitons, as well as an area ("unstable solitons") where the unstable antisymmetric 
soliton spontaneously transforms itself into a single-peak mode trapped in either circle, and shedding off about a half 
of the initial norm. In terms of the evolution of the respective density profile, \ip (x,y,t)\ 2 , the latter outcome of the 
evolution is very similar to that shown in Fig. [3] 

Note that only parts of the solution branches with the negative slope, dN/d/u, < 0, which meet the VK criterion, 
are stable in Fig. [?ta) (see the curves pertaining to a = 1, 1.5, 2, and 3). The stability and instability intervals for 
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FIG. 4: (Color online) (a) An example of the transformation of a weakly unstable symmetric soliton into a robust breather, 
with N — 22, R = 4.3n and a = 1.5. (b) The evolution of largest values of the local density of the left ("L") and right ("R") 
peaks, which features a weak spontaneous breaking of the symmetry between the peaks. 




FIG. 5: (Color online) The spontaneous transformation of an unstable symmetric double-peak soliton into a single-peak one, 
at N = 22, R = 4tt and a = 1.2. 
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a 



FIG. 6: (Color online) The stability diagram for various localized modes developing from stationary symmetric solitons at a 
fixed value of the norm, N = 22, in the plane of the radius of the circles (a) and normalized separation between their centers 
(R/tt). "Unstable solitons" are those which spontaneously transform into the single-peak modes, as shown in Fig. [5] The 
dashed line represents L — 0, i.e., R = 2a. Along this line, the two circles are tangent to each other, and below it the circles 
partly overlap. As shown in the next section, only single-peak solitons may be stable at L < 0. 




H N 

(a) (b) 

FIG. 7: (Color online) The same as in Fig. [2] (for the fixed separation between the circles, L = 37r), but for families of 
antisymmetric stationary solitons. 

these three values of a are reported in Table 1. Further, similar to the situation for the symmetric solitons displayed 
in Fig. H[b), the E(N) curves for the antisymmetric solutions displayed in Fig. EJb) confirm the natural expectation 
that stable portions of the solution families have lower energy than their unstable counterparts, for a = 1 and 1.5. 

C. Overlapping and touching circles (L < 0) 

As shown in Figs. [6] and [51 stable double-peak solitons, symmetric or antisymmetric ones, cannot be supported by 
the pair of circles with a small separation between them. Further numerical analysis demonstrates that such solitons 
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FIG. 8: (Color online) The same as in Fig. [6](with the fixed norm N = 22), but for antisymmetric solitons and modes generated 
by the development of their instability. 

TABLE I: Stability intervals, in terms of the total norm (N) and chemical potential (fi), for families of stationary antisymmetric 
solitons presented in Fig. for fixed L — 3tt. 



a 


N 


M 




1.0 


21.8347 < N < 22.4702 


-1.5 < fj, < 


-0.68 


1.5 


21.9289 < N < 23.0877 


-1.2 < n < 


-0.42 


2.0 


22.1718 < N < 23.2130 


-0.8 < (i < 


-0.3 


3.0 


22.9914 < N < 23.3356 


-0.5 < y, < - 


-0.278 



cannot be found either for touching (L — 0) or overlapping (L < 0) circles. Instead, in this case it is possible to find 
symmetric single-peak solitons, centered at the midpoint of the configuration. Strongly asymmetric solitons, trapped 
near the center of either circle, while the other one is left almost empty, exist in this case too, and they may be stable. 
These modes are considered below. 



1. Touching circles (L = 0) 

In the case of L = 0, the two circles touch each other at the single point, featuring a "figure of eight" [see insets in 
Fig. |9ja)]. As said above, in this case symmetric stationary modes are represented by single-peak solitons centered at 
the touch point, see examples of the soliton profiles in Fig. [Hta). These solitons are unstable, spontaneously jumping 
into one of the circles, as seen in Fig. [2Jb) . This instability can be easily understood, as the symmetric soliton is 
actually located at the position of an unstable equilibrium between two attractive pseudopotential wells corresponding 
to the circles. After the jump, the single-peak soliton keeps hopping in the circle, periodically hitting its border and 
bouncing back. 

Strongly asymmetric single-peak solitons trapped near the center of either circle can be found too at L — 0. An 
example of such a stable nearly isotropic soliton in shown in Fig. 1101 The numerical data demonstrate that the 
ratio of local densities at the center of the circle, and at the midpoint where the two circles touch each other, is 
[(f) (1.5, 0) /((> (0, 0)] 2 ~ 25 for this soliton. 

Figure [TT1 displays N(fj,) and E(N) dependences for families of such single-peak asymmetric modes, obtained at 
different values of radius a. It is seen that stable subfamilies are found for 1 < a < 3. The stability of these 
solution families exactly obeys the VK criterion, with unstable modes featuring decay into radiation (not shown here) . 
Extremely close to the critical value (|2"4"|) . the stationary single-peak solitons develop a weak oscillatory instability 
and turn into breathers featuring intrinsic regular vibrations with a small amplitude (not shown here in detail). For 
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(a) (b) 

FIG. 9: (Color online) (a) Profiles of unstable symmetric single-peak solitons centered around the touch point of the set of two 
circles with zero separation, L = (see the insets). The top and bottom panels pertain to different radii of the circles, a — 2 
and a — 4. In both cases, the norm of the soliton is TV = 15.9. (b) The ensuing hopping motion of the single-peak breather 
inside one circle, shown by means of the density contour plots in the x-cross section. 




n 



-4 
3 
■2 



FIG. 10: (Color online) An example of the stable single-peak soliton, trapped in the right circle of the "figure-of-eight" structure 
(with L = 0), is shown by means of contour plots of the density, (<j> (x,y)) 2 , for fi = —1, N = 11.4798, and a = 1.5. 



instance, at a = 1.5, the robust breathers are observed in an interval of 11.668 < N < 11.690. 



D. Partly overlapping circles (L < 0) 

Negative L (—2a < L < 0) implies that the two circles partly overlap [see the top left inset in Fig. [T27a)]. merging 
into a single circle at L = —2a. Only single-peak solitons were found in this case. First, we continue the numerical 
analysis of the symmetric solitons centered at the midpoint, x — y — 0. In a narrow interval corresponding to 
slightly overlapping circles, —0.1a < L < 0, the symmetric solitons are unstable in essentially the same fashion as 
demonstrated above in Fig. [9l i.e., they jump from the unstable equilibrium position into either circle, performing 
the shuttle motion in it. 

As the degree of the overlap between the circles increases, in the interval of —0.33a < L < —0.1a the unstable 
symmetric soliton simply decays into radiation (not shown here in detail). This outcome of the evolution may be 
interpreted as a consequence of the fact that the total norm of the soliton is too low in this case, not allowing it 
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FIG. 11: (Color online) The same as in Figs. [2] but for the family of stationary single-peak solitons nested in either circle 
of the "figure-of-eight" configuration (the one with L — 0). As before, the red color designates stable solutions, while black 
branches represent modes which are unstable through the decay into radiation. 



to maintain its integrity. Further, keeping to increase the overlap degree, in the interval of —0.33a < L < 0.48 we 
observe that the gradual enhancement of the nonlinearity around the midpoint does not yet stabilize the symmetric 
soliton, but restores its integrity. In this case, the soliton again spontaneously leaps to the left or right circle, which 
is followed by its shuttle motion inside the circle. 

The symmetric solitons centered around the midpoint may be stable at L < —0.48a. For instance, at L = —0.50a, 
they are stable in the interval of iV m i n = 11.2523 < N < 11.2665 = N maX) which corresponds to —0.55 < fi < —0.47, 
cf. Table 1. Similarly, at L = -0.55a the stability interval is 11.2375 < N < 11.2637, or -0.59 < n < -0.47. The 
stability interval for the symmetric solitons expands with the increase of \L\ towards L — —2a. 

The situation when the symmetric solitons may be stable is illustrated by Fig. [T2] The right inset to panel [T27a) 
demonstrates that, at fi = —0.6 (N = 11.2592), the N(/j) curve splits into two branches, which then merge at 
/i = —1.65 (N — 11.5295). The resulting bifurcation loop resembles those found in the ID and 2D double-core models 
with the cubic-quintic nonlinearity [TH, [yj . The splitting means that two soliton solutions with different values of the 
norm can be found at given /i (for instance, the pair of solitons corresponding to points "c" and "d" in Fig. [T2")) . In 
the splitting region, the top branch in panel IT2T a) represents the symmetric soliton placed at the midpoint, while the 
lower branch represents a soliton shifted to left or right. In fact, the splitting and recombination of the two solution 
branches are explicit examples of the direct and inverse symmetry-breaking bifurcations. 

Note that the solitons of types "a" and "d", shown in Fig. IT2T b). are centered exactly at the midpoint, while the 
solitons of types "b" and "c" show a small shift off the midpoint. Of course, together with the shifted soliton, it is 
possible to find its mirror-image counterpart, shifted in the opposite direction. 

The simulations demonstrate that broad symmetric solitons, corresponding to large values of TV [cf. Eq. (|2"0"|) ]. suffer 
decay at n > —0.47 (i.e., at N > 11.2361), as indicated by the dashed black line in Fig. IT2T a). Stable symmetric 
solitons where found in the interval of —0.6 < [i < —0.47 (in terms of the norm, it is 11.2654 < N < 11.2361), 
which corresponds to the short solid red segment in Fig. IT2T a). At /i < —0.6, the solitons marked by "b" and "c", 
which belong to the family indicated by the dashed-dotted magenta curve, spontaneously transform themselves into 
breathers which feature regular small-amplitude oscillations. Further, solitons pertaining to the dotted blue curve in 
Fig. fT2Ta) (for instance, the soliton corresponding to point "d" in Fig. [T2")) undergo a spontaneous transformation 
into breathers which feature strong irregular intrinsic oscillations, but remain robust localized modes. 

The splitting of the N(fi) curve occurs in a small region of the plane of (a, L/a), as shown in Fig. Q21 As concerns 
the (narrow) stability area of the symmetric solitons, it is outlined by Table 2, for the overlapping parameter fixed 
with respect to the radius of the circles, L = —0.5720a (the same relation between L and a as in Fig. [T2")) . At values of 
N exceeding the largest value indicated in Table 2, the symmetric solitons spontaneously turn into robust single-peak 
breathers. 

Finally, for strongly overlapping circles, with —2a ^ L < —0.8, the results, including the stability region for the 
symmetric solitons (only the symmetric ones exist in this case) are nearly the same as reported in Ref. |36j for the 
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(a) (b) 

FIG. 12: (Color online) (a) The N((j,) curves for stationary single-peak solitons in the configuration with L = —0.5720a and 
a = 1.1. In part (b), subplots labeled "a", "b", "c" and "d" display cross sections along x (the magenta dashed lines) and along 
y (the blue solid lines) of the single-peak solitons corresponding to points "a", "b", "c" and "d" in panel (a), i.e., respectively, 
for (j, = -0.56, N = 11.2544; fx = -0.6, N = 11.2592; fi = -0.8, N = 11.2606; and fi = -0.8, N = 11.3421. 








Ua 



-1.5 



_2" ■ ■ ■ 
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a 

FIG. 13: (Color online) The region in the plane of (a,L/a), at fixed /i = —0.8, where curves N(fi) for the family of single-peak 
solitons, supported by the overlapping circles, feature the splitting [see Fig. fTST a - )]. 

single circle, which corresponds to L — —2a. 

V. SUMMARY 

In this work, we have introduced the 2D model based on the modulation function confining the self-focusing cubic 
nonlinearity, i.e., the corresponding nonlinear potential, to the two identical circles, which may be separated or 
overlapped. We aimed to study symmetric, antisymmetric, and asymmetric stationary 2D solitons in this setting, as 
well as breathers into which unstable solitons may be spontaneously transformed. Results of the analysis, obtained 
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TABLE II: Stability intervals for families of stationary symmetric solitons centered around the midpoint of the configuration 
with partly overlapping circles, L — —0.5720a, and several values of radius a of the overlapping circles. 



a 


N 


A* 




1.1 


11.2361 < N < 11.2620 


-0.47 < a* < 


-0.588 


1.5 


11.3449 < N < 11.2633 


-0.33 < p, < 


-0.47 


1.8 


11.2025 < N < 11.2611 


-0.17 < p, < 


-0.25 


a > 2 


no stability intervals 







by means of numerical methods and the VA (variational approximation), can be summarized as follows below. 

Well-separated circles support stable symmetric and antisymmetric solitons, as long as the interaction between the 
density peaks trapped in the two circles is negligible. The increase of the interaction strength (i.e., the decrease of 
separation L between the circles, relative to their radii a) leads to the destabilization of the symmetric and antisym- 
metric solitons. As shown in Fig. [6l the symmetric solitons pass two SSB transitions. At first, they are transformed 
into robust breathers featuring small-amplitude irregular intrinsic vibrations, with a spontaneously emerging small 
asymmetry between the two peaks (Fig. @|. As the interaction between the two circles grows still stronger, the weakly 
asymmetric double-peak breathers lose their stability and are replaced by single-peak modes residing in one circle, 
while the other one remains nearly empty. On the other hand, Fig. [8] shows that the antisymmetric solitons feature 
the single SSB transition, by a direct jump into the single-peak mode. The VA, based on the superposition of two 
Gaussians centered in the two circles, describes the symmetric solitons reasonably well, while it fails to accurately 
approximate antisymmetric ones. 

The situation is different in the case of touching (L = 0) and overlapping (L < 0) circles. In this case, only single- 
peak solitons are found - both asymmetric modes, trapped in one circle (Fig. I10j) . a part of which are stable (Fig. 
Hip , and symmetric solitons centered around the midpoint (x = y = 0). For the weak overlapping, < — L < 0.1a, 
the symmetric single-peak soliton is situated at the unstable-equilibrium position. As a result, it spontaneously leaps 
into either circle and then performs shuttle motion therein (Fig. |9]). For a stronger overlap, 0.1a < — L < 0.33a, the 
symmetric soliton decays into radiation. At a still larger degree of the overlap, 0.48a < — L < 0.33a, the enhancement 
of the local nonlinearity around the midpoint results in the reappearance of the SSB regime, in the form of the 
spontaneous leap followed by the shuttle motion. A region of the stability of the symmetric single-peak solitons 
appears at —L > 0.48a, gradu ally expanding, with the growth of \L\, into the stability region for the soliton trapped 
in the single nonlinear circle [36| in the limit when the two circles merge into the single one (— L — > 2a). In the 
case of the moderately strong overlap, a noteworthy manifestation of the SSB was found in the form of the paired 
symmetry-breaking and symmetry- restoring bifurcations, which give rise to the loop in Fig. 112( a). 

This work may be naturally extended in other directions. In particular, a related dynamical problem is Josephson 
oscillations in bosonic junctions formed by double- well potentials, cf. Refs. Q, (4^-|4a|. The present model calls for 
the study of Josephson oscillations in the double-well nonlinear pseudopotentials. Because the setting has no linear 
limit, the Josephson frequency is expected to strongly depend on the norm of the condensate. The models considered 
in Ref. [33| and here suggest that the corresponding nonlinear Josephson junctions may be experimentally realized 
in photonic media, in addition to BEC. 

The nonlinear pseudopotentials may be also naturally generalized for two-component models, which give rise to 
vector solitons [35]]. Accordingly, it may be relevant to study SSB effects in double- well pseudopotentials trapping 
two-component mixtures. 
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